% This Matlab program will conduct the bootstrap estimations necessary to generate figure 4;

clear
tic

load ForBootstraps

reps = 10000;

stcnty = double( st_dataset.stcnty );
n = size(stcnty,1);
S = sparse( [1:292308]', stcnty, ones(n,1) );  % (n x m) '

tram = double(st_dataset.tr_am);
availcredit = double(st_dataset.availcredit);
limit = double(st_dataset.limit);
pctblack = double(st_dataset.pctblack);
util_dol = double(st_dataset.util_dol);
bal = double(st_dataset.bal);
sc_black = double(st_dataset.sc_black);
unusedlim = double(st_dataset.unusedlim);

demographics = [ st_dataset.age st_dataset.gt_hs_male st_dataset.eq_hs_male st_dataset.gt_hs_female st_dataset.eq_hs_female st_dataset.nonmarried_male st_dataset.nonmarried_female st_dataset.widowed_male st_dataset.widowed_female st_dataset.divorced_male st_dataset.divorced_female st_dataset.foreignborn ];
income = [ st_dataset.income_10 st_dataset.income_15 st_dataset.income_20 st_dataset.income_25 st_dataset.income_30 st_dataset.income_35 st_dataset.income_40 st_dataset.income_45 st_dataset.income_50 st_dataset.income_75 st_dataset.income_100 st_dataset.income_125 st_dataset.income_150 st_dataset.income_200 ];
emp = [ st_dataset.publicassistance st_dataset.employment ];
housing = [ st_dataset.vacant st_dataset.ownocc st_dataset.withmort st_dataset.medianrent_b st_dataset.mhv ];

clear st_dataset;

n = size(availcredit,1);
[sort_fips,fips_index] = sort(stcnty);
chg_fips = zeros(n,1);
chg_fips(1) = 1;
chg_fips(2:n) = (sort_fips(2:n)~=sort_fips(1:n-1));
fips_grp = cumsum(chg_fips,1);
fips_grp(fips_index) = fips_grp;
fips_dummies = sparse(dummyvar(fips_grp));
%fips_dummies(:,1) = [];

% Baseline with interaction - AVAILCREDIT
y = availcredit;
x = [pctblack sc_black tram];

[n,k] = size(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)
fits = x*beta;
residuals = y - fits;
penalty = - x(:,1:2)*beta(1:2);
mynodes = linspace(0,1,100);
my_penalty = kreg(x(:,1),penalty,0.05,mynodes);

% Now, let's bootstrap the standard errors;
node_mat = zeros(100,reps);
parfor count=1:reps;
		if (floor(count/10)*10 == count)
			fprintf('%d\n',count);
		end;
    epsilon = residuals(ceil(n.*rand(n,1)));
    newy = fits + epsilon;
    newbeta = xpxi*(x'*newy);  % '
    
    penalty = - x(:,1:2)*newbeta(1:2);
    node_mat(:,count) = kreg(x(:,1),penalty,0.05,mynodes);
end;
sort_penalty = sort(node_mat,2);
ci_low = sort_penalty(:,floor(reps*0.025));
ci_high = sort_penalty(:,ceil(reps*0.975));

figure(1);
clf;
set(gcf, 'PaperUnits','inches','PaperPosition',[0 0 6.5 9]);

subplot(3,2,1);
fill([mynodes fliplr(mynodes)],[ci_high' fliplr(ci_low')],[0.85 0.85 0.85],'LineStyle','none');
title('(a) AVAILCREDIT Replication');
xlabel('PCT\_BLACK');
ylabel('Race Penalty ($000s)');
axis([0 1 0 12]);
hold on;
plot(mynodes,my_penalty,'k',mynodes,zeros(100,1),'k-');
hold off;
set(gca,'position',[0.1 0.72 .35 .25])
holda = [mynodes' my_penalty ci_high ci_low];

% Baseline with interaction - LIMIT
y = limit;
x = [pctblack sc_black util_dol tram];

[n,k] = size(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)
fits = x*beta;
residuals = y - fits;
penalty = - x(:,1:2)*beta(1:2);
mynodes = linspace(0,1,100);
my_penalty = kreg(x(:,1),penalty,0.05,mynodes);

% Now, let's bootstrap the standard errors;
node_mat = zeros(100,reps);
parfor count=1:reps;
		if (floor(count/10)*10 == count)
			fprintf('%d\n',count);
		end;
    epsilon = residuals(ceil(n.*rand(n,1)));
    newy = fits + epsilon;
    newbeta = xpxi*(x'*newy);  % '
    
    penalty = - x(:,1:2)*newbeta(1:2);
    node_mat(:,count) = kreg(x(:,1),penalty,0.05,mynodes);
end;
sort_penalty = sort(node_mat,2);
ci_low = sort_penalty(:,floor(reps*0.025));
ci_high = sort_penalty(:,ceil(reps*0.975));

subplot(3,2,2);
fill([mynodes fliplr(mynodes)],[ci_high' fliplr(ci_low')],[0.85 0.85 0.85],'LineStyle','none');
title('(b) LIMIT Replication');
xlabel('PCT\_BLACK');
ylabel('Race Penalty ($000s)');
axis([0 1 0 12]);
hold on;
plot(mynodes,my_penalty,'k',mynodes,zeros(100,1),'k-');
hold off;
set(gca,'position',[.57 .72 .35 .25])
holdb = [mynodes' my_penalty ci_high ci_low];


% Baseline with interaction and income without distortion - AVAILCREDIT
y = unusedlim;
x = [pctblack sc_black tram income];

[n,k] = size(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)
fits = x*beta;
residuals = y - fits;
penalty = - x(:,1:2)*beta(1:2);
mynodes = linspace(0,1,100);
my_penalty = kreg(x(:,1),penalty,0.05,mynodes);

% Now, let's bootstrap the standard errors;
node_mat = zeros(100,reps);
parfor count=1:reps;
		if (floor(count/10)*10 == count)
			fprintf('%d\n',count);
		end;
    epsilon = residuals(ceil(n.*rand(n,1)));
    newy = fits + epsilon;
    newbeta = xpxi*(x'*newy);  % '
    
    penalty = - x(:,1:2)*newbeta(1:2);
    node_mat(:,count) = kreg(x(:,1),penalty,0.05,mynodes);
end;
sort_penalty = sort(node_mat,2);
ci_low = sort_penalty(:,floor(reps*0.025));
ci_high = sort_penalty(:,ceil(reps*0.975));

subplot(3,2,3);
fill([mynodes fliplr(mynodes)],[ci_high' fliplr(ci_low')],[0.85 0.85 0.85],'LineStyle','none');
title('(c) AVAILCREDIT Without Distortion Effect Plus Income');
xlabel('PCT\_BLACK');
ylabel('Race Penalty ($000s)');
axis([0 1 -6 2]);
hold on;
plot(mynodes,my_penalty,'k',mynodes,zeros(100,1),'k-');
hold off;
set(gca,'position',[0.1 0.39 .35 .25])
holdc = [mynodes' my_penalty ci_high ci_low];


% Baseline with interaction and income without distortion - LIMIT
y = limit;
x = [pctblack sc_black bal tram income];

[n,k] = size(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)
fits = x*beta;
residuals = y - fits;
penalty = - x(:,1:2)*beta(1:2);
mynodes = linspace(0,1,100);
my_penalty = kreg(x(:,1),penalty,0.05,mynodes);

% Now, let's bootstrap the standard errors;
node_mat = zeros(100,reps);
parfor count=1:reps;
		if (floor(count/10)*10 == count)
			fprintf('%d\n',count);
		end;
    epsilon = residuals(ceil(n.*rand(n,1)));
    newy = fits + epsilon;
    newbeta = xpxi*(x'*newy);  % '
    
    penalty = - x(:,1:2)*newbeta(1:2);
    node_mat(:,count) = kreg(x(:,1),penalty,0.05,mynodes);
end;
sort_penalty = sort(node_mat,2);
ci_low = sort_penalty(:,floor(reps*0.025));
ci_high = sort_penalty(:,ceil(reps*0.975));

subplot(3,2,4);
fill([mynodes fliplr(mynodes)],[ci_high' fliplr(ci_low')],[0.85 0.85 0.85],'LineStyle','none');
title('D) AVAILCREDIT Without Distortion Effect Plus Income');
xlabel('PCT\_BLACK');
ylabel('Race Penalty ($000s)');
axis([0 1 -6 2]);
hold on;
plot(mynodes,my_penalty,'k',mynodes,zeros(100,1),'k-');
hold off;
set(gca,'position',[.57 .39 .35 .25])
holdd = [mynodes' my_penalty ci_high ci_low];


% Baseline with interaction, income, and other controls without distortion - AVAILCREDIT
y = unusedlim;
x = [pctblack sc_black tram income demographics emp housing];

[n,k] = size(x);
x = double(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)
fits = x*beta;
residuals = y - fits;
penalty = - x(:,1:2)*beta(1:2);
mynodes = linspace(0,1,100);
my_penalty = kreg(x(:,1),penalty,0.05,mynodes);

% Now, let's bootstrap the standard errors;
node_mat = zeros(100,reps);
parfor count=1:reps;
		if (floor(count/10)*10 == count)
			fprintf('%d\n',count);
		end;
    epsilon = residuals(ceil(n.*rand(n,1)));
    newy = fits + epsilon;
    newbeta = xpxi*(x'*newy);  % '
    
    penalty = - x(:,1:2)*newbeta(1:2);
    node_mat(:,count) = kreg(x(:,1),penalty,0.05,mynodes);
end;
sort_penalty = sort(node_mat,2);
ci_low = sort_penalty(:,floor(reps*0.025));
ci_high = sort_penalty(:,ceil(reps*0.975));

subplot(3,2,5);
fill([mynodes fliplr(mynodes)],[ci_high' fliplr(ci_low')],[0.85 0.85 0.85],'LineStyle','none');
title('E) AVAILCREDIT Without Distortion Effects Plus Full Controls');
xlabel('PCT\_BLACK');
ylabel('Race Penalty ($000s)');
axis([0 1 -6 2]);
hold on;
plot(mynodes,my_penalty,'k',mynodes,zeros(100,1),'k-');
hold off;
set(gca,'position',[0.1 0.06 .35 .25])
holde = [mynodes' my_penalty ci_high ci_low];


% Baseline with interaction, income, and other controls without distortion - LIMIT
y = limit;
x = [pctblack sc_black bal tram income demographics emp housing];

[n,k] = size(x);
x = double(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)
fits = x*beta;
residuals = y - fits;
penalty = - x(:,1:2)*beta(1:2);
mynodes = linspace(0,1,100);
my_penalty = kreg(x(:,1),penalty,0.05,mynodes);

% Now, let's bootstrap the standard errors;
node_mat = zeros(100,reps);
parfor count=1:reps;
		if (floor(count/10)*10 == count)
			fprintf('%d\n',count);
		end;
    epsilon = residuals(ceil(n.*rand(n,1)));
    newy = fits + epsilon;
    newbeta = xpxi*(x'*newy);  % '
    
    penalty = - x(:,1:2)*newbeta(1:2);
    node_mat(:,count) = kreg(x(:,1),penalty,0.05,mynodes);
end;
sort_penalty = sort(node_mat,2);
ci_low = sort_penalty(:,floor(reps*0.025));
ci_high = sort_penalty(:,ceil(reps*0.975));

subplot(3,2,6);
fill([mynodes fliplr(mynodes)],[ci_high' fliplr(ci_low')],[0.85 0.85 0.85],'LineStyle','none');
title('F) LIMIT Without Distortion Effects, Plus Full Controls');
xlabel('PCT\_BLACK');
ylabel('Race Penalty ($000s)');
axis([0 1 -6 2]);
hold on;
plot(mynodes,my_penalty,'k',mynodes,zeros(100,1),'k-');
hold off;
set(gca,'position',[.57 .06 .35 .25])
holdf = [mynodes' my_penalty ci_high ci_low];

% This generages the graphic file that becomes figure 4;
print -dpng NewBootstrapGraph.png -r100

save('BootReps.mat','holda','holdb','holdc','holdd','holde','holdf');

toc

